DGE analysis of mPFC
Libraries required
Load salmon object
Differeitial analysis using limma
All samples
design <- model.matrix(~ 0 + Group, data = salmon@phenoData)
colnames(design) <- gsub(pattern = "Group", replacement = "", x = colnames(design))
y <- DGEList(counts = salmon@gene.counts)
keep <- filterByExpr(y, design, min.count = 20)
y <- y[keep, ]
dds <- calcNormFactors(y)
v <- voom(dds, design = design)
contrast.matrix <- makeContrasts(MSUS - CTRL, levels = design)
fit <- lmFit(v)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
dea <- as.data.frame(topTable(fit2, coef = 1, number = Inf))
dea <- data.frame(Genes = rownames(dea), dea, stringsAsFactors = F)
dea_sig <- dea[abs(dea$logFC) >= 1 & dea$adj.P.Val <= 0.05, ]Removing L001
design_l1 <- model.matrix(~ 0 + Group, data = salmon@phenoData[-1,])
colnames(design_l1) <- gsub(pattern = "Group", replacement = "", x = colnames(design_l1))
y_l1 <- DGEList(counts = salmon@gene.counts[,-1])
keep_l1 <- filterByExpr(y_l1, design_l1, min.count = 20)
y_l1 <- y_l1[keep_l1, ]
dds_l1 <- calcNormFactors(y_l1)
v_l1 <- voom(dds_l1, design = design_l1)
contrast.matrix_l1 <- makeContrasts(MSUS - CTRL, levels = design_l1)
fit_l1 <- lmFit(v_l1)
fit2_l1 <- contrasts.fit(fit_l1, contrast.matrix_l1)
fit2_l1 <- eBayes(fit2_l1)
dea_l1 <- as.data.frame(topTable(fit2_l1, coef = 1, number = Inf))
dea_l1 <- data.frame(Genes = rownames(dea_l1), dea_l1, stringsAsFactors = F)
dea_sig_l1 <- dea_l1[abs(dea_l1$logFC) >= 1 & dea_l1$adj.P.Val <= 0.05, ]cpm Counts data and pData
Results
dea.list <- list(`CTRL vs MSUS` = as.DEA(dea),
`CTRL vs MSUS (no L001)` = as.DEA(dea_l1))
normCounts <- v$E
voomEList <- v
dea.limma <- list(`CTRL vs MSUS` = dea)Save RData files
save(dea.list,
file = "./output/dea_mPFC_lit.DEA.RData", compress = T,
compression_level = 3
)
save(normCounts,
file = "./output/normCounts_voom_mPFC_lit.RData", compress = T,
compression_level = 3
)
save(voomEList,
file = "./output/voom_EList_mPFC_lit.RData", compress = T,
compression_level = 3
)
save(dea,
file = "./output/limma_mPFC_lit.RData", compress = T,
compression_level = 3
)
writexl::write_xlsx(x = dea, path = "output/dea_results.xlsx", col_names = T, format_headers = T)MA plots
All samples
res <- dea
res$threshold <- as.factor(res$adj.P.Val < 0.1)
p1 <- ggplot(data = res, aes(
x = res$AveExpr,
y = res$logFC,
colour = threshold
)) +
geom_point(alpha = 0.5, size = 1.8) +
geom_hline(aes(yintercept = 0), colour = "blue", size = 1) +
ylim(c(
-ceiling(max(abs(res$logFC))),
ceiling(max(abs(res$logFC)))
)) +
ggtitle("MSUS vs CTRL") +
labs(subtitle = "loess fit") +
xlab("Mean expression") +
ylab("Log2 Fold Change") +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
plot.subtitle = element_text(size = 16, hjust = 0.5, face = "italic", color = "blue"),
axis.title = element_text(face = "bold", size = 15),
axis.text = element_text(face = "bold", size = 12),
legend.title = element_text(face = "bold", size = 15),
legend.text = element_text(size = 14),
panel.border = element_rect(colour = "black", fill=NA, size=2)
) +
scale_color_manual(name = "p.adjusted < 0.1", values=brewer.pal(n = 3, name = "Set1")[1:2])+ NULL
stat_smooth(se = FALSE, method = "loess", color = "red", formula = y ~ x, size = 1)## geom_smooth: se = FALSE, na.rm = FALSE, orientation = NA
## stat_smooth: method = loess, formula = y ~ x, se = FALSE, n = 80, fullrange = FALSE, level = 0.95, na.rm = FALSE, orientation = NA, method.args = list(), span = 0.75
## position_identity
Removing L0001
res_l1 <- dea_l1
res_l1$threshold <- as.factor(res_l1$adj.P.Val < 0.1)
p2 <- ggplot(data = res_l1, aes(
x = res_l1$AveExpr,
y = res_l1$logFC,
colour = threshold
)) +
geom_point(alpha = 0.8, size = 3) +
geom_hline(aes(yintercept = 0), colour = "blue", size = 1) +
ylim(c(
-ceiling(max(abs(res_l1$logFC))),
ceiling(max(abs(res_l1$logFC)))
)) +
ggtitle("MSUS vs CTRL (no L1)") +
labs(subtitle = "loess fit") +
xlab("Mean expression") +
ylab("Log2 Fold Change") +
theme_bw() +
theme(
plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
plot.subtitle = element_text(size = 16, hjust = 0.5, face = "italic", color = "blue"),
axis.title = element_text(face = "bold", size = 15),
axis.text = element_text(face = "bold", size = 12),
legend.title = element_text(face = "bold", size = 15),
legend.text = element_text(size = 14),
panel.border = element_rect(colour = "black", fill=NA, size=2)
) +
scale_color_manual(name = "p.adjusted < 0.1", values=brewer.pal(n = 3, name = "Set1")[1:2])+ NULL
stat_smooth(se = FALSE, method = "loess", color = "red", formula = y ~ x, size = 1)## geom_smooth: se = FALSE, na.rm = FALSE, orientation = NA
## stat_smooth: method = loess, formula = y ~ x, se = FALSE, n = 80, fullrange = FALSE, level = 0.95, na.rm = FALSE, orientation = NA, method.args = list(), span = 0.75
## position_identity
## Warning: Use of `res$AveExpr` is discouraged. Use `AveExpr` instead.
## Warning: Use of `res$logFC` is discouraged. Use `logFC` instead.
## Warning: Use of `res_l1$AveExpr` is discouraged. Use `AveExpr` instead.
## Warning: Use of `res_l1$logFC` is discouraged. Use `logFC` instead.
Heatmap of most variable genes
All samples
xvar <- apply((data$cpm + 1), 1, var)
genes500 <- head(sort(xvar, decreasing = TRUE), n = 500)
x500 <- data$cpm[rownames(data$cpm) %in% names(genes500), ]
pheatmap(x500, scale = "row", show_rownames = F, col = viridis(100), show_colnames = F, annotation_col = data$pData[,2:3], main = "Clustering of samples based on top 500 variable genes")Removing L001
xvar <- apply((data$cpm[,-1] + 1), 1, var)
genes500 <- head(sort(xvar, decreasing = TRUE), n = 500)
x500 <- data$cpm[rownames(data$cpm) %in% names(genes500), ]
pheatmap(x500, scale = "row", show_rownames = F, col = viridis(100), show_colnames = F, annotation_col = data$pData[,2:3], main = "Clustering of samples based on top 500 variable genes (no L1)")PCA
All samples
pca <- prcomp(t(data$cpm))
# ggplotly(
pc1 <- autoplot(pca,
data = data$pData, colour = "Group",
frame = TRUE, frame.type = "norm", size = 5
) +
ggtitle("") +
theme_bw() +
scale_color_manual(values=brewer.pal(n = 3, name = "Set1")[1:2]) +
theme(axis.text.x = element_text(size = 20, colour = "black"),
axis.text.y = element_text(size = 20, colour = "black"),
axis.title.x = element_text(size = 20, colour = "black"),
axis.title.y = element_text(size = 20, colour = "black"),
plot.title = element_text(size = 25, face = "bold", colour = "black"),
legend.text=element_text(size=20, colour = "black"),
legend.title=element_text(size=20, face = "bold"),
panel.border = element_rect(colour = "black", fill=NA, size=2))
pc2 <- autoplot(pca,
data = data$pData, colour = "Lane",
frame = TRUE, frame.type = "norm", size = 5
) +
ggtitle("") +
theme_bw()+
scale_color_manual(values=brewer.pal(n = 3, name = "Set1")[1:2]) +
theme(axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20),
plot.title = element_text(size = 25, face = "bold"),
legend.text=element_text(size=20),
legend.title=element_text(size=20, face = "bold"),
panel.border = element_rect(colour = "black", fill=NA, size=2))
pca1 <- prcomp(t(data$cpm[,-1]))
# ggplotly(
pc3 <- autoplot(pca1,
data = data$pData[-1,], colour = "Group",
frame = TRUE, frame.type = "norm", size = 5
) +
ggtitle("") +
theme_bw() +
scale_color_manual(values=brewer.pal(n = 3, name = "Set1")[1:2]) +
theme(axis.text.x = element_text(size = 20, colour = "black"),
axis.text.y = element_text(size = 20, colour = "black"),
axis.title.x = element_text(size = 20, colour = "black"),
axis.title.y = element_text(size = 20, colour = "black"),
plot.title = element_text(size = 25, face = "bold", colour = "black"),
legend.text=element_text(size=20, colour = "black"),
legend.title=element_text(size=20, face = "bold"),
panel.border = element_rect(colour = "black", fill=NA, size=2))## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
All samples
##
## Attaching package: 'plotly'
## The following object is masked from 'package:IRanges':
##
## slice
## The following object is masked from 'package:S4Vectors':
##
## rename
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning: Didn't find a colorbar to modify.
Remove L1
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning: Didn't find a colorbar to modify.
References
## Warning in citation(pkg_name): no date field in DESCRIPTION file of package
## 'plgINS'
SessionInfo
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 3.6.2 (2019-12-12)
os Ubuntu 16.04.7 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate de_DE.UTF-8
ctype de_DE.UTF-8
tz Europe/Zurich
date 2020-11-10
─ Packages ───────────────────────────────────────────────────────────────────
package * version date lib
annotate 1.64.0 2019-10-29 [1]
AnnotationDbi 1.48.0 2019-10-29 [1]
assertthat 0.2.1 2019-03-21 [1]
autoplotly * 0.1.2 2018-04-21 [1]
backports 1.1.8 2020-06-17 [1]
base64enc 0.1-3 2015-07-28 [1]
Biobase * 2.46.0 2019-10-29 [1]
BiocGenerics * 0.32.0 2019-10-29 [1]
BiocManager 1.30.10 2019-11-16 [1]
BiocParallel * 1.20.1 2019-12-21 [1]
Biostrings 2.54.0 2019-10-29 [1]
bit 4.0.4 2020-08-04 [1]
bit64 4.0.2 2020-07-30 [1]
bitops 1.0-6 2013-08-17 [1]
blob 1.2.1 2020-01-20 [1]
bookdown 0.20 2020-06-23 [1]
callr 3.4.3 2020-03-28 [1]
checkmate 2.0.0 2020-02-06 [1]
cli 2.0.2 2020-02-28 [1]
cluster 2.1.0 2019-06-19 [1]
colorspace 1.4-1 2019-03-18 [1]
crayon 1.3.4 2017-09-16 [1]
crosstalk 1.1.0.1 2020-03-13 [1]
data.table 1.13.0 2020-07-24 [1]
DBI 1.1.0 2019-12-15 [1]
DelayedArray * 0.12.3 2020-04-09 [1]
desc 1.2.0 2018-05-01 [1]
DESeq2 1.26.0 2019-10-29 [1]
devtools 2.3.1 2020-07-21 [1]
digest 0.6.25 2020-02-23 [1]
dplyr * 1.0.1 2020-07-31 [1]
edgeR * 3.28.1 2020-02-26 [1]
ellipsis 0.3.1 2020-05-15 [1]
evaluate 0.14 2019-05-28 [1]
fansi 0.4.1 2020-01-08 [1]
farver 2.0.3 2020-01-16 [1]
foreign 0.8-76 2020-03-03 [1]
Formula 1.2-3 2018-05-03 [1]
fs 1.5.0 2020-07-31 [1]
genefilter * 1.68.0 2019-10-29 [1]
geneplotter 1.64.0 2019-10-29 [1]
generics 0.0.2 2018-11-29 [1]
GenomeInfoDb * 1.22.1 2020-03-27 [1]
GenomeInfoDbData 1.2.2 2019-11-18 [1]
GenomicRanges * 1.38.0 2019-10-29 [1]
GEOquery 2.54.1 2019-11-18 [1]
ggfortify 0.4.10 2020-07-27 [1]
ggplot2 * 3.3.2 2020-06-19 [1]
ggplotify * 0.0.5 2020-03-12 [1]
glue 1.4.1 2020-05-13 [1]
gridExtra 2.3 2017-09-09 [1]
gridGraphics 0.5-0 2020-02-25 [1]
gtable 0.3.0 2019-03-25 [1]
Hmisc 4.4-1 2020-08-10 [1]
hms 0.5.3 2020-01-08 [1]
htmlTable 2.0.1 2020-07-05 [1]
htmltools 0.5.0 2020-06-16 [1]
htmlwidgets 1.5.1 2019-10-08 [1]
httr 1.4.2 2020-07-20 [1]
insight 0.9.0 2020-07-20 [1]
IRanges * 2.20.2 2020-01-13 [1]
jpeg 0.1-8.1 2019-10-24 [1]
jsonlite 1.7.0 2020-06-25 [1]
knitr 1.29 2020-06-23 [1]
labeling 0.3 2014-08-23 [1]
lattice 0.20-41 2020-04-02 [1]
latticeExtra 0.6-29 2019-12-19 [1]
lazyeval 0.2.2 2019-03-15 [1]
lifecycle 0.2.0 2020-03-06 [1]
limma * 3.42.2 2020-02-03 [1]
locfit 1.5-9.4 2020-03-25 [1]
magrittr 1.5 2014-11-22 [1]
Matrix 1.2-18 2019-11-27 [1]
matrixStats * 0.56.0 2020-03-13 [1]
memoise 1.1.0.9000 2020-05-06 [1]
mgcv * 1.8-31 2019-11-09 [1]
munsell 0.5.0 2018-06-12 [1]
nlme * 3.1-148 2020-05-24 [1]
nnet 7.3-14 2020-04-26 [1]
patchwork * 1.0.1 2020-06-22 [1]
pheatmap * 1.0.12 2019-01-04 [1]
pillar 1.4.6 2020-07-10 [1]
pkgbuild 1.1.0 2020-07-13 [1]
pkgconfig 2.0.3 2019-09-22 [1]
pkgload 1.1.0 2020-05-29 [1]
plgINS * 0.1.5 2020-07-10 [1]
plotly * 4.9.2.1 2020-04-04 [1]
plyr 1.8.6 2020-03-03 [1]
png 0.1-7 2013-12-03 [1]
prettyunits 1.1.1 2020-01-24 [1]
processx 3.4.3 2020-07-05 [1]
ps 1.3.3 2020-05-08 [1]
purrr 0.3.4 2020-04-17 [1]
R6 2.4.1 2019-11-12 [1]
RColorBrewer * 1.1-2 2014-12-07 [1]
Rcpp 1.0.5 2020-07-06 [1]
RCurl 1.98-1.2 2020-04-18 [1]
readr 1.3.1 2018-12-21 [1]
remotes 2.2.0 2020-07-21 [1]
report 0.1.0 2020-03-19 [1]
reshape2 * 1.4.4 2020-04-09 [1]
rlang 0.4.7 2020-07-09 [1]
rmarkdown 2.3 2020-06-18 [1]
rmdformats 0.4.0 2020-06-07 [1]
rpart 4.1-15 2019-04-12 [1]
rprojroot 1.3-2 2018-01-03 [1]
RSQLite 2.1.4 2019-12-04 [1]
rstudioapi 0.11 2020-02-07 [1]
rvcheck 0.1.8 2020-03-01 [1]
S4Vectors * 0.24.4 2020-04-09 [1]
scales 1.1.1 2020-05-11 [1]
sessioninfo 1.1.1 2018-11-05 [1]
SRAdb 1.48.2 2019-12-24 [1]
stringi 1.4.6 2020-02-17 [1]
stringr 1.4.0 2019-02-10 [1]
SummarizedExperiment * 1.16.1 2019-12-19 [1]
survival 3.2-3 2020-06-13 [1]
sva * 3.34.0 2019-10-29 [1]
testthat 2.3.2 2020-03-02 [1]
tibble 3.0.3 2020-07-10 [1]
tidyr 1.1.0 2020-05-20 [1]
tidyselect 1.1.0 2020-05-11 [1]
usethis 1.6.1 2020-04-29 [1]
vctrs 0.3.2 2020-07-15 [1]
viridis * 0.5.1 2018-03-29 [1]
viridisLite * 0.3.0 2018-02-01 [1]
withr 2.2.0 2020-04-20 [1]
writexl 1.3 2020-05-05 [1]
xfun 0.16 2020-07-24 [1]
XML 3.99-0.3 2020-01-20 [1]
xml2 1.3.2 2020-04-23 [1]
xtable 1.8-4 2019-04-21 [1]
XVector 0.26.0 2019-10-29 [1]
yaml 2.2.1 2020-02-01 [1]
zlibbioc 1.32.0 2019-10-29 [1]
source
Bioconductor
Bioconductor
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
Bioconductor
CRAN (R 3.6.1)
Bioconductor
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
Bioconductor
Bioconductor
CRAN (R 3.6.1)
Bioconductor
Bioconductor
Bioconductor
Bioconductor
Github (sinhrks/ggfortify@3f4020d)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
Github (r-lib/memoise@4aefd9f)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
local
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
Github (easystats/report@dcdd283)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
Github (juba/rmdformats@94cd7a3)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.2)
Bioconductor
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
CRAN (R 3.6.1)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.2)
CRAN (R 3.6.1)
Bioconductor
CRAN (R 3.6.2)
Bioconductor
[1] /home/ubuntu/R/x86_64-pc-linux-gnu-library/3.6
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library